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! Abstract 
^ . 

' We introduce a one-dimensional toy model of globular clusters. The model is 

a version of the well-known gravitational sheets system, where we take addition- 
ally into account mass and energy loss by evaporation of stars at the boundaries. 
' Numerical integration by the "exact" event-driven dynamics is performed, for ini- 

• tial uniform density and Gaussian random velocities. Two distinct quasi-stationary 

asymptotic regimes are attained, depending on the initial energy of the system. We 
O ' guess the forms of the density and velocity profiles which fit numerical data ex- 

-4— > ' tremely well and allow to perform an independent calculation of the self-consistent 

, gravitational potential. Some power-laws for the asymptotic number of stars and 

^ I for the collision times are suggested. 

X ; PACS numbers: 05.45.-a; 05.45.Pq; 98.10.+Z; 98.20.Jp. 

P.: 

1 Introduction 

Globular clusters are gravitationally bound concentrations of large numbers of stars, 
spherically distributed in space. They orbit around a galaxy spending most of the time in 
the galactic halo 1^. The most important elements governing globular clusters structure 
are two-body relaxation and truncation due to tidal forces. Different dynamical models 
considering these specific phenomena, have been investigated both analytically |jl|, |], |^, Q 
and numerically [|, |, 0. For an extended discussion, see and Refs. therein. 
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Dynamical evolution causes stars to escape as an effect of the gravitational interaction 
with the nearby galaxies. This evaporation process drives the cluster towards a configura- 
tion with a high density core and the velocity dispersion of stars in the bulk can increase 
without limit. This phenomenon is known as gravothermal catastrophe and its study goes 
back to Antonov and to Lynden-Bell and Wood |]I0[ . 

Referring to the pioneering work of Chandrasekhar it is possible to calculate the 
perturbations induced by stellar encounters on star motion. This is done by means of a 
diffusion model, which allows to derive a quantitative description of changes of star veloc- 
ities in terms of single encounters. Considering weak encounters, i.e. solving the diffusion 
model in the Fokker-Planck approximation. King found the following expression for 
the velocity profile 



where Vc is a cutoff velocity of the stars, a is the one-dimensional velocity dispersion and 
A a normalization constant. King's models have been shown to be in agreement with 
the observed brightness surface profiles of globular clusters |Q. Further developments 
of King's model p[ consider the same functional form (truncated Gaussian) applied to 
the gravitational energy, and hence propose a general form for the distribution function 
/(x, v). We will not consider here these extensions. 

In this paper we discuss a simplified one-dimensional A^-body model which reproduces 
King's distribution. We let the particles, all of equal mass m, interact through the one- 
dimensional gravitational potential V = 27rG'm^|x|, where G is the gravitational constant. 
Bearing in mind the comparison with globular clusters, we imitate the effect of galactic 
tidal forces by introducing a finite cut-off in positions. Thus, the evaporation of stars 
from the system is the only "dissipative" effect we consider. This is enough to drive the 
system towards an asymptotic non stationary regime, which we analyze in detail, and 
which reveals striking similarity with King's model. 

The main difference between the simplified one-dimensional model considered here 
and a more realistic three-dimensional one, is the lack of any singularity of the potential 
at the origin. The presence of a finite lower bound for the ID potential makes less 
energy available to support the evaporation process as the system cools down. In the 
3D case an infinite amount of energy can indeed be extracted from the singular pair- wise 
gravitational interaction, which is the main origin of the gravothermal catastrophe. This 
is the reason why the model we discuss cannot reproduce the core collapse corresponding 
to the gravothermal catastrophe. 

In the next Section we present the model and the results concerning velocity distri- 
bution and density profiles. Section ^ is devoted to the discussion of power-laws for the 
number of particles in the cluster. Finally, in Section ^ we draw some conclusions. 
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2 One-dimensional model 



Let us consider a one- dimensional classical Newtonian self-gravitating system of N parti- 



cles with equal mass m, with Hamiltonian [|ri 

1 ^ 

H = -m^Vi + 2nGm^ ^\xi-Xj\ , (2) 

i j<i 

where Xi is the position and Vi the velocity of i-th star. G is the universal gravitational 
constant. We choose in the following m = 1 and 27iG = 1. This system has been recently 



the subject of intensive investigations |jT2[. Particle accelerations are constant inbetween 
two collisions and are proportional to the net difference of particles respectively on the 
right and on the left. When a collision occurs particles cross each other, or, equivalently, 
collisions are elastic. This particle approach is known to corresponds in the continuum 
limit (A^ oo) to the Vlasov-Poisson equations for the distribution function f{x,v) 

^ + ^^_^^ = 
dt dx dx dv 



dx"^ 



AnGm J f(x,v)dv , (3) 



where V{x) is the self-consistent gravitational potential. 
We add the following features to the model 

• Particles are confined in a box of size L, i.e. Xi G [—L/2, L/2] 

• The effect of tidal forces induced by the parent galaxy is imitated by requiring 
that each time a particle reaches the boundary of the box with a finite velocity, 
it drops out of the system, which therefore experiences a mass and energy loss 
("evaporation"). This last feature implies that system (|^) is solved with absorbing 
boundary conditions. 

The numerical implementation is based on an "event-driven" scheme, first introduced 
in plasma physics [|13[, which is adapted to the present follows. The algorithm 



looks for the particles which collide the first and for the time when the event occurs, tcou- 
Then it computes the first "evaporation time", tevap, and makes the system evolve until 
the minimum, tmin, between the collision and the evaporation time is reached. Once the 
system experiences evaporation, the total mass is reduced and the escaping particle stops 
interacting with the residual bulk. By rescaling the position and the velocity of each 
particle, i.e. introducing a local dissipation, we maintain the position of the center of 
mass fixed and its velocity to zero. This means we simply translate velocity and position 
of the remaining particles to keep the system centered in position and momentum space. 
A particle can escape from the system as a result of this rescaling: this possibility is taken 
into account even if it is has a low probability. 
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Evaporation is a singular event, which, in fact, marks the transition between two self- 
gravitating system having a different number of particles and energy. We remark that the 
integration scheme is "exact". Time t elapsed from the initial configuration is obtained 
by summing all values of tmin up to the last event. 

In all numerical experiments, N{0) particles are initially uniformly distributed in the 
box and the initial velocities are Gaussian i.i.d random variables, with the temperature 
To given by twice the average kinetic energy. 

In Fig. |l] we show the time evolution of the number of particles N{t) which remain 
inside the box up to time t. After an abrupt decrease of N(t), which strongly depends 
on the initial condition, the system reaches a state where rare evaporations are present, 
making N{t) decrease much slowlier. Our numerical experiments show that such a quasi- 
stationary state lives indefinitely, although we cannot exclude that, finally, N{t) relaxes 
to an asymptotic value Nas- As the best approximation for this value, we take the one 
the system reaches in the longest computer runs. 

Given L, for large enough Tq the system approaches a state characterized by a single 
cluster, which adapts itself to the size of the box. In Fig. |^ we show the phase-space 
portrait for a system in its late stage evolution, when the most energetic particles have 
dropped out from the box and the system has relaxed to an asymptotic plateau, as the 
ones reported in Fig. |l|. The particles are almost uniformly distributed within a bounded 
region of the phase plane. Both this fact and the shape of the contour suggest a possible 



connection with the so called water-bag (WB) distribution [jll] , |14|. This is a station- 
ary solution of the Vlasov-Poisson system f^^{x,v), which is constant in a simply 
connected domain Q of the phase-plane and strictly zero outside. 

Adopting the notation of Ref. fll]], we call respectively Xg and Vg the maximum position 
and velocity of the WB. The potential V{x) for such a distribution is then implicitly 
specified by the following integral equation 
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X 



Vix) 



ei-(6-C)^ 'dC. (4) 



The maximal energy of the water-bag, e, is such that, if we express f{x, v) in terms of the 
energy u = + V{x), f{x, v) = F{u) = if u > e. The energy e is related to Xg and 
f s by e = fs/2, V{xs) = e and the zero energy level is fixed by requiring that ^(0) = 0. 
The density profile, p^^(x), is expressed as function of the potential 



Since the distribution function f{x,v) is constant over Q, it follows immediately 

fix,v)dv = 2cv+{x) (6) 

-v+{x) 

where v^{x) represents the profile of the upper branch of the WB contour, which we 
assumed to be symmetric, and c = SN^^/ {16^/26^^"^). Using Eq. (|^), this implies 

v+{x) = v/2(e - V) . (7) 
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To compute the velocity contour v~^{x), we need to know V{x), which we do by solving 
Eq. (Q) by an adaptive recursive Newton-Cotes 8 panel rule with tolerance 10~^. This 
velocity contour is drawn in Fig. ^ and, as predicted, encloses all phase points. In deriving 
V{x) we have taken the value Vg from the cluster phase plot; this is the only phenomeno- 
logical input in this calculation and the agreement with the data has to be considered 
quite satisfactory. 

Moreover, we can compare the theoretically derived potential V{x) with the one com- 
puted directly from the asymptotic positions 

Vixi) ='^\xi- Xj\ . (8) 
j 

We need, of course, to perform a vertical shift to fix the zero in the origin. The result 
of formula @ is reported in Fig. ^ together with the theoretically derived potential. 
The agreement is really good. We are thus led to conclude that our asymptotic state 
is well described by a water-bag. However, this latter is a stationary solution of the 
self-gravitating ID system, while in our simulations we continue to observe some particle 
evaporations even at very long times. This is why we have called our asymptotic state 
quasi- stationary and its description in terms of a water-bag distribution can only be 
approximate. 

An alternative treatment of the asymptotic state is based on King's formula (|l]). In this 
case one does not try to reproduce the full distribution function, but just its projections 
along the x and v axis: p(x) = / fdv and fo{v) = J fdx. Following the standard 
derivation of the equilibrium isothermal distribution |T^, we are led to introduce an 
analytical ansatz for the density profile 



p{x) = A{cosh~'^{Bx) - cosh~^{B^)) for |a;| < 



p(x) = for |x| > ^ 



(9) 



L 

2 ' 



T In 

where the normalization A is fixed by J_^^^ p{x)dx = Nas- For the velocities we take 
King's distribution, specified in Eq. (|1|). Assuming p(x) as in Eq. (j^) we can derive a 
close analytical expression for the potential. For a one- dimensional system, the following 
relation holds in the continuum limit 

V{x) = J'jy~x\p{y)dy . (10) 
Inserting Eq. (^) into Eq. ([T0|) and performing the integral, we get 



V{x) = Vo 



— tanh(5-) - ix^ + — ) cosh"^ fi-) + — ln( 

B ^ 2' ^ a' ^ 2' B^ cosh(5f) 



'111 
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with 

^° " |tanh(5f) -Lcosh-2(5f) ' ^^^^ 

is quadratic for small x. To verify the reliability of our guess we re-analyze the data 
previously discussed in connection with the water-bag distribution. In Fig. ^ we plot the 
potential calculated numerically from Eq. (^) together with a one-parameter fit, using 
Eq. (0). Again the agreement with the data is very good and even apparently superior 
to the one obtained using the WB picture. This is simply due to the fact that here we 
perform a one parameter fit, while, in the previous discussion, Vs was arbitrarily deduced 
from the phase space analysis. As a cross-check, we introduce in Eq. (|^) the coefficient B 
determined from the fit of the potential. The resulting density profile is plotted in Fig. |^ 
and it agrees with the normalized histogram of particles position. Finally, an histogram 
of the velocity is represented in Fig. ^. The reverse-cup shape due to the cut-off of the 
tails is evident. The solid line in Fig. ^ is a numerical fit which uses the expression of 
Eq. (|1|) with Vc and a as free parameters. 

As a side remark we observe that, coherently with the observed form of the potential, 
each particle oscillates almost harmonically inside the box. This can be seen by looking 
at the asymptotic orbit of a single particle (Fig. |^). The slight diffusion of the orbit is 
the signature of the interaction with the other particles, which induces a weak chaoticity. 

A further aspect that we have tested is the dependence of the dynamics on the initial 
temperature, Tq. Indeed, for small values of Tq, the system shows a pronounced collapse, 
which leads to a massive central core, as it is clearly displayed in the main plot of Fig. ^. 
This phase-space distribution significantly differs from the one in Fig. ^ and cannot be 
represented by a water-bag. The histograms of positions and velocities are computed 
and plotted in the right and left inserts, respectively. Both the density and the velocity 
profiles are very well reproduced by a numerical fit based on our ansatz (y) and on King's 
distribution (|I]). 

In conclusion, the ansatz we have introduced shows a good agreement with numerical 
data for all values of Tq we have simulated, while the water-bag distribution fails to 
reproduce the velocity and density profiles at very low initial temperature. However, in 
the high temperature range, the water-bag treatment is superior, because it leads to an 
accurate description of the full distribution f{x,v). 



3 Scaling laws 

In this Section we discuss some numerically found scaling laws which do not have presently 
a theoretical justification but which are an important signature of the presence of a finite 
box. 

We follow the system until it reaches the asymptotic state. Being Ncoii the number of 
collisions, we define the "average collision time", r = tt^- 

In the main plot of Fig we represent Nas as a function of r. Each point refers to a 
different value of the initial temperature To, varying form 0.2 to 7, while Ncou is maintained 
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constant for each realization. The initial temperature controls the rate of evaporation at a 
very early stage of the evolution. Larger values of Tq, produce higher mass loss, inducing 
the system to relax to a quasi- stationary state characterized by less residual particles Nas 
(see Fig. |l]). This process has consequences at the dynamical level, determining a larger 
mean free path and consequently a larger value of r. The curve in Fig ^ is consistent 
with this qualitative picture, showing a power-law decay with exponents a = —0.4 which 
is valid over more than two decades. In the insert of Fig. |^, we plot r vs. Tq. 

We have also checked the dependence on N of the collision time tcoii, defined as an 
average of all the collision times corresponding to each fixed value of A^, i.e. inbetween two 
successive evaporations. In Fig. |10| we plot the results of different numerical experiments 
where we vary the initial temperature Tq. The dependence of tcoii on N is again a power- 
law with exponent b = —2.5. Note that b ~ 1/a, hence Fig. ^ can be thought as a 



macroscopic averaged image of the microscopic properties shown in Fig. 10 



4 Conclusions 

We have introduced a one- dimensional toy model of globular clusters with an emphasis on 
the evaporation process. With this in mind, we have discussed the effect of introducing a 
finite size box in a classical one dimensional self-gravitating medium. The dynamics of the 
system has been investigated for a special class of initial conditions. We pointed out the 
appearance of two distinct, non stationary, asymptotic regimes which occur depending 
on the temperature of the initial realization. For small values of Tq, similarities with 
the isothermal solution are found while for larger temperatures the density and velocity 
profiles are well reproduced also assuming a water-bag distribution. 

We propose a form of the density profile, with a cut-off in the tails, which fits well 
numerical data in all the explored regimes, allowing to derive a close analytical expression 
of the gravitational potential. Moreover, a King-like velocity profile is shown to be in 
good agreement with numerical data. The asymptotic truncated profiles are thus a direct 
consequence of the evaporation from the finite box. 

Finally we have also given numerical evidence of some scaling laws, which remain to 
be theoretically explained, but which are strongly related to the escaping process. 

In the future, we plan to extend this study to the system of concentric spherical mass 
shells by introducing an external absorbing boundary in the configuration space as 
done here. 
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Figure 1: Plots of vs. time for increasing initial temperatures, with A^(0) = 400. 
Temperature and time are expressed in arbitrary units. 




Figure 2: Phase-space plot for a system of A^(0) = 1500 particles after 30 x 10® collisions. 
Here Tq = 0.4, L = 0.0015 and Nas = 886. The quantities x and v are expressed in 
arbitrary units. The full line which contains all the points is the theoretical prediction 

m 
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Figure 3: Gravitational potential calculated numerically using Eq. (squares) and 
analytically by solving the integral equation (D with Vg = 0.88 (full line). Only the 
region of positive x is drawn, in arbitrary units. 
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Figure 4: Gravitational potential calculated numerically using Eq. (H) (squares) and by 
a numerical fit which uses Eq. ( [TT| ) (full line), where B = 731.2 is the only free parameter. 
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Figure 5: Normalized histogram of positions as derived from the phase-space plot in 
Fig. 1^. The solid line is Eq. @, with B = 731.2 as previously. Position x is expressed in 
arbitrary units. 
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Figure 6: Normalized histogram of velocities as derived from the phase-space plot in 

Fig. |. The solid line is the fit obtained using Eq. (|ID with cr = 0.87, Vc = 0.98. Velocity 
V is expressed in arbitrary units. 
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Figure 7: Asymptotic orbit of a single particle for iV(0) = 600, L = 0.0035. Positions 
and velocities are expressed in arbitrary units. 




Figure 8: Phase-space plot for a system of A^(0) = 1500 particles after 30 x 10® collisions. 
Here Tq = 0.02, L = 0.0015 and Nas = 1167. Left insert: normalized histogram of 
positions. The solid line is a fit which uses Eq. (§) where B = 6845.7. Right insert: 
normalized histogram of velocities. The solid line is the fit obtained using Eq. ([^) with 
a = 0.3, Vc = 1.9. All quantities are expressed in arbitrary units. 
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Figure 9: Nas vs. r in log- log scale. Here r is defined as the ratio 77^- Each point 
refers to a different Tq while Ncoii is fixed. The solid line represents a power-law fit with 
the slope a = —0.4. In the upper right corner insert r vs. Tq is represented in a log-log 
scale. The quantities Tq and r are expressed in arbitrary units. 
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Figure 10: tcoii vs. N in log-log scale for A^(0) = 600 L = 0.0035. Different symbols refer 
to different initial values of the temperature Tq. tcou is expressed in arbitrary units. 
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